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Abstract 

We  present  a  (suboptimal)  filtering  algorithm  for  tracking  a  highly  maneuvering  target  in  a 
cluttered  environment  using  multiple  sensors.  The  filtering  algorithm  is  developed  by  applying 
the  basic  Interacting  Multiple  Model  (IMM)  approach  and  the  Probabilistic  Data  Association 
(PDA)  technique  to  a  two  sensor  (radar  and  infrared,  for  instance)  problem  for  state  estimation 
for  the  target.  A  detection  fusion  approach  is  followed  where  the  raw  sensor  measurements  are 
passed  to  a  fusion  node  and  fed  directly  to  the  target  tracker.  A  multisensor  probabilistic  data 
association  filter  is  developed  for  parallel  sensor  processing  for  target  tracking  under  clutter.  A 
past  approach  using  parallel  sensor  processing  has  ignored  certain  data  association  probabilities 
leading  to  an  erroneous  derivation.  Another  existing  approach  applies  only  to  non-maneuvering 
targets.  The  algorithm  is  illustrated  via  a  highly  maneuvering  target  tracking  simulation  ex¬ 
ample  where  two  sensors,  a  radar  and  an  infrared  sensor,  are  used.  Compared  with  an  existing 
IMMPDA  filtering  algorithm  with  sequential  sensor  processing,  the  proposed  algorithm  achieves 
significant  improvement  in  the  accuracy  of  track  estimation. 
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1  Introduction 


We  consider  the  problem  of  tracking  a  single  maneuvering  target  in  clutter.  This  class  of  problem 
has  received  considerable  attention  in  the  literature  [1,  6,  7,  8,  9,  12].  The  switching  multiple  model 
approach  has  been  found  to  be  quite  effective  in  modeling  highly  maneuvering  targets  [1,  2,  3,  4, 
5,  6,  9].  In  this  approach  various  ’’modes”  of  target  motion  are  represented  by  distinct  kinematic 
models,  and  in  a  Bayesian  framework,  the  target  maneuvers  are  modeled  by  switching  among  these 
models  controlled  by  a  Markov  chain.  In  the  presence  of  clutter,  the  measurements  at  the  sensors 
may  not  all  have  originated  from  the  target-of-interest.  In  this  case  one  has  to  solve  the  problem 
of  data  association.  An  effective  approach  in  a  Bayesian  framework  is  that  of  probabilistic  data 
association  (PDA)  [2,  12]. 

To  accommodate  the  fact  that  the  target  can  be  highly  maneuvering,  we  will  follow  a  switching 
multiple  model  formulation  as  in  [1,  2,  3,  4,  5,  6]  and  references  therein.  It  is  assumed  that  a 
track  has  been  formed  (initiated)  and  our  objective  is  that  of  track  maintenance.  In  [1]  such  a 
problem  has  been  considered  using  multiple  sensors,  PDA,  and  switching  multiple  models.  The 
optimal  solution  (in  the  minimum  mean-square  error  sense)  to  target  state  estimation  given  sensor 
measurements  and  absence  of  clutter,  requires  exponentially  increasing  (with  time)  computational 
complexity;  therefore,  one  has  to  resort  to  suboptimal  approximations.  For  the  switching  multiple 
model  approach,  the  interacting  multiple  model  (IMM)  algorithm  of  [4]  has  been  found  to  offer  a 
good  compromise  between  the  computational  and  storage  requirements  and  estimation  accuracy 
[3].  In  the  presence  of  clutter,  one  has  to  account  for  measurements  of  uncertain  origin  (target 
or  clutter?).  Here  too,  in  a  Bayesian  framework,  one  has  to  resort  to  approximations  to  reduce 
the  computational  complexity,  resulting  in  the  PDA  filter  [1,  6,  8,  9,  12].  In  [1]  the  IMM  algo¬ 
rithm  has  been  combined  with  the  PDA  filter  in  a  multiple  sensor  scenario  to  propose  a  combined 
IMM/MSPDAF  (interacting  multiple  model/  multisensor  probabilistic  data  association  filter)  al¬ 
gorithm.  The  multisensor  approach  of  [1]  falls  in  the  category  of  detection  fusion  where  the  raw 
measurements  from  all  sensors  are  passed  to  a  fusion  node  to  be  processed  simultaneously  [9].  Ref 
[1]  uses  sequential  updating  of  the  state  estimates  with  measurements  (i.e.  updating  the  state 
estimates  sequentially  with  measurements  from  different  sensors).  This  results  in  computational 
savings  but  this  approach  is  not  necessarily  the  best.  The  other  option  is  that  of  parallel  updating 
(i.e.  updating  the  state  estimates  with  all  the  measurements  at  the  same  time  as  if  they  were  from 
a  single  sensor).  For  linear  systems,  the  two  updating  methods  are  algebraically  equivalent  but 
the  parallel  updating  is  computationally  more  expensive  [6].  Ref  [9]  uses  parallel  updating  but  has 
some  errors:  during  data  association,  all  measurements  at  the  same  time  from  different  sensors  are 
assumed  to  be  either  from  clutter  of  from  the  target.  The  possibility  that  a  measurement  from 
sensor  1  may  be  from  target  while  the  measurement  from  sensor  2  may  be  clutter-induced  (and 
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vice-versa)  in  implicitly  not  allowed  in  [9]  -  this  is  clearly  incorrect.  Ref  [13]  allows  for  such  dis¬ 
tinctions  (hypotheses),  however,  it  is  limited  to  non-maneuvering  targets.  In  this  paper,  we  extend 
the  multisensor  approach  of  [13]  to  maneuvering  targets. 

The  paper  is  organized  as  follows.  Section  2  presents  the  problem  formulation.  Section  3 
describes  the  proposed  IMM/MSPDAF  algorithm  with  parallel  detection  fusion.  Simulation  results 
using  the  proposed  algorithm  for  a  realistic  problem  are  given  in  Section  4.  Finally,  Section  5 
presents  a  discussion  of  the  results  and  some  conclusions. 

2  Problem  Formulation 

Assume  that  the  target  dynamics  can  be  modeled  by  one  of  n  hypothesis  models.  The  model  set 
is  denoted  as  A4n  :=  {1,  ...,n}  and  there  are  tot^  q  sensors  from  which  5,  or  fewer  (if  probability 
of  target  detection  is  less  than  one)  or  more  (due  to  clutter),  measurement  vectors  are  generated 
at  a  time.  The  event  that  model  j  is  in  effect  during  the  sampling  period  is  denoted  by 

M^.  For  the  y-th  hypothesized  model  (mode),  the  state  dynamics  and  measurements,  respectively, 
modeled  as 

Xk  =  Fk-i^k-i  +  (1) 

and 

zi=h^’''{xk)  +  wi:^  for  (2) 

where  Xk  is  the  system  state  at  tk  and  of  dimension  Ux,  is  the  (true)  measurement  vector  (i.e.,  due 
to  the  target)  from  sensor  I  at  tk  and  of  dimension  Uzi,  Fl_^  and  are  the  system  matrices  when 
model  j  is  in  effect  over  the  sampling  period  (tfc_i,tA;],  and  is  the  nonlinear  transformation  of 
Xk  to  z\.  {I  —  1, ...,  9)  for  model  j.  Henceforth,  time  tk  will  be  denoted  hy  k.  A  first-order  linearized 
version  of  (2)  is  given  by 

Zk  =  Hi'^Xk  +  w];^  for  l  =  l,...,q  (3) 

where  is  the  Jacobian  matrix  of  evaluated  at  some  value  of  the  estimate  of  state  Xk-  The 
process  noise  and  the  measurement  noise  are  mutually  uncorrelated  zero-mean  white 
Gaussian  processes  with  covariance  matrices  Qk_i  and  Iti’\  respectively.  At  the  initial  time  to, 
the  initial  conditions  for  the  system  state  under  each  model  j  are  assumed  to  be  Gaussian  random 
variables  with  the  known  mean  ^  and  the  known  covariance  Pj-  The  probability  of  model  j  at 
^o^  Mo  ~  F{Mq},  is  also  assumed  to  be  known.  The  switching  from  model  to  model 

is  governed  by  a  finite-state  stationary  Markov  chain  with  known  transition  probabilities  pij  = 
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The  following  notations  and  definitions  are  used  regarding  the  measurements  at  sensor  1.  Note 
that,  in  general,  at  any  time  k,  some  measurements  may  be  due  to  clutter  and  some  due  to  the 
target,  i.e.  there  can  be  more  than  a  single  measmement  at  time  k  at  sensor  1.  The  measurement 
set  (not  yet  validated)  generated  by  sensor  I  at  time  k  is  denoted  as 


zi 


1(1)  1(2) 


} 


(4) 


where  mi  is  the  number  of  measurements  generated  by  sensor  I  at  time  k.  Variable  (z  = 
is  the  ith  measurement  within  the  set.  The  validated  set  of  measurements  of  sensor  I  at 
time  k  will  be  denoted  by  Yl^  containing  (<  m/)  measurement  vectors.  The  cumulative  set  of 
validated  measurements  firom  sensor  I  up  to  time  k  is  denoted  as 


(5) 


The  cumulative  set  of  validated  measurements  from  all  sensors  up  to  time  k  is  denoted  as 

Zk  {^fc(l)^  ^fc(2)^  ^*(9)  j 

where  q  is  the  number  of  sensors. 

Our  goal  is  to  find  the  state  estimate 

Xk\k-=^E{xk\Z>^}  (7) 

and  the  associated  error  covariance  matrix 


Ek\k  • —  \Z  }  (8) 

where  x'j^  denotes  the  transpose  of  Xk- 

3  IMM/MSPDAF  Algorithm  for  Parallel  Detection  Fusion 

We  now  modify  the  IMM/PDA  algorithms  of  [1]  and  [11]  to  derive  the  proposed  IMM/MSPDAF 
with  parallel  detection  fusion  system.  We  confine  our  attention  to  the  case  of  2  sensors;  however, 
the  algorithm  can  be  easily  adapted  to  the  case  of  arbitrary  q  sensors.  We  will  only  briefly  outline 
the  basic  steps  in  “one  cycle”  (i.e.,  processing  needed  to  update  for  a  new  set  of  measurements)  of 
the  IMM/MSPDA  filter. 

Assumed  available:  Given  the  state  estimate  :=  E{xk-i\Ml_.y,  the  associated 

covariance  conditional  mode  probability  :=  at  time  /i;  —  1  for 

each  mode  j  E  M.n- 
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Step  1.  Interaction  —  mixing  of  the  estimate  from  the  previous  time  (ij  G  M.n)  • 
predicted  mode  probability: 

i 

mixing  probability: 

:=  P[MU\Mi  =  Py4-i/4“-  (10) 

mixed  estimate: 

•=  E{x,.,\Mi,Z'‘-')  =  (11) 

i 

covariance  of  the  mixed  estimate: 

^k-l\k-l  •“  ^{[^k-l  — 

t 

Step  2.  Predicted  state  and  measurements  for  sensors  1  and  2  (ij  G  Mn)  : 
state  prediction: 

4|.-l  Z‘-'}  -fl'-iCilt-r  (13) 

state  prediction  error  covariance: 


The  mode-conditioned  predicted  measurement  for  sensor  I  is 


Using  the  linearized  version  (3),  the  covariance  of  the  mode-conditioned  residual 


(14) 

(15) 


is  given  by  (assume  g=2,  the  case  of  2  sensors) 

Si'  :=  Z‘-‘)  ;=  Hi'%-l»f  +  ^4'“.  (16) 

where  is  the  first  order  derivative  (Jacobian  matrix)  of  h^'^{.)  evaluated  at  the  state  prediction 
^k\k-i  (®®®  (10))-  Note  that  (16)  and  (17)  assume  that  z^''^  originates  from  the  target.  The  results 
(16)  and  (17)  do  not  depend  upon  the  actual  measurements. 
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As  mentioned  earlier,  since  our  approach  to  the  problem  deals  with  multiple  simultaneous 
measurements  [13,  14]  arising  from  two  separate  sensors  that  are  tracking  a  single  target  through 
a  common  surveillance  region,  a  method  for  fusion  of  multiple  measurements  has  to  be  devised.  In 
order  to  do  this,  now  the  combined  covariance  5^  of  the  mode-conditioned  residual  obtained  from 
(16)  and  (17)  also  needs  to  be  considered  as  follows: 


Step  3.  Measurement  validation  for  sensors  1  and  2  (ij  6  Mn)- 
There  is  uncertainty  regarding  the  measurements’  origins.  Therefore,  we  perform  validation  for 
each  target  separately.  One  sets  up  a  validation  gate  for  sensor  I  centered  at  the  mode-conditioned 
predicted  measurement,  Let  (|A|  =  det(A.)) 

Then  measurement  (i=:l,2,...,mi)  is  validated  if  and  only  if 

(19) 

where  7  is  an  appropriate  threshold.  The  volume  of  the  validation  region  with  the  threshold  7  is 

vt  ■■=  (20) 

where  n^i  is  the  dimension  of  the  measurement  and  Cn^^  is  the  volume  of  the  unit  hypersphere  of 
this  dimension  (ci  =  2,  C2  =  tt,  C3  =  47r/3,  etc.).  Choice  of  7  is  discussed  in  more  detail  in  ([6],  Sec. 
2.3.2).  After  performing  the  validation  for  each  target  separately,  we  deal  with  all  the  validated 
data  for  the  measurement  fusion. 


Step  4'  State  estimation  with  validated  measurement  from  sensors  1  and  2  (Vj  G 
Mn)  : 

Prom  among  all  the  raw  measurements  from  sensor  I  at  time  A:,  i.e.  ;=  ...,2^”*'^}, 

define  the  set  of  validated  measurement  for  sensor  I  at  time  k  as 


ii  ~  {vF'.yf* . %'*'’} 


where  mi  is  total  number  of  validated  measurement  for  sensor  I  at  time  k  and 


._  Aii) 
Vk  •- 


where  1<  li  <  I2  <  ...  <  /m,  <  fni  when  mi  ^0.  Define  the  association  events  (hypotheses)  9’^^  as 
follows  (here  we  follow  [?]) 
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•  none  of  the  measurements  in  or  is  target  originated. 

•  6^^:  only  in  Y^  is  a  target  measurement,  all  other  measurements  in  Y^  or  Y^  are  clutter, 

o  =  0,  h  =  1, 

•  6^'^:  only  in  Y^  is  a  target  measurement,  all  other  measurements  in  Y^  or  Y^  are  clutter, 
a  =  1, ...,  mi,  b  =  0. 

•  0^’*’:  and  in  Y/^  and  Y^,  respectively,  are  target  measurements,  all  other  measure¬ 
ments  are  clutter,  a  =  1,  6  =  1,  ...,m2. 

Therefore,  there  are  a  total  of  mim2+mi+m2+l  possible  association  hypotheses,  each  of  which 
has  an  association  probability.  Define  the  mode-conditioned  association  event  probabilities  as 


Exploiting  the  diffuse  model  for  clutter  in  [1,  6],  it  turns  out  that 

Ojfi,0  ^  )(j-~ ^02-^02)  _  ri  L  n 

P/c  ^  ^  (ylNmj  (Y2\m2  ^ ^  ®  ^ 


(23) 


\j,0,b  _  Q 


=  C 


Pni(l-PDjPG2)N 


m2 


,  0  =  0,  6=1,  ...,m2 


o4'‘ 


(24) 


fh2 


3^'“’*  =  c 


N 


,  0  =  1,  ...,mi,  6  =  0 
]  Pd,Pd2 


0  =  1, ..., mi,  6  =  1,..., m2 


where  Pd^  and  P02  are  the  detection  probabilities  that  the  sensors  1  and  2  detect  the  target, 
respectively,  Pq^  and  Pci  are  probabilities  the  target  is  in  the  validation  region  observed  from 
sensors  1  and  2,  respectively,  C  is  a  normalization  constant  such  that  ~  1  '^3 


N[x;y,P]  :=  |27rP|  ^'^^exp 
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Define  the  mode-conditioned  innovations  as 


= 


Onzixi 

0n22xi 

Oflilxl 


,  a  =  0,  b  =  0 


,  a  =  0,  6  =  1, 


■  j/Ma) 

^  ^  a  =  l,...,mi,  b  =  0 

0fl22xl 

r  ^,l(a)  - 

*4’“’*=  ,  a  =  b=  1,. 


,  a  =  b=  l,...,m2. 


The  likelihood  function  for  each  mode  m  is 


mi  7712 

4  ■■=v\yiX\Miz'‘-']  =  y.Y.p 


a=0 6=0 


where 


a  =  0,  b  =  0 


X  N  ] ,  a  =  0,  b  =  l,...,m2 

=  1  °2L  fci  ^27) 

X  N  ]  ,  a=l . mi,(.  =  0 

.  ^  ^  ■  “  =  1 . 

Using  (from  (13))  and  its  covariance  (^o™  (l'^));  computes  the  partial  update 

and  its  covariance  P^f.  according  to  the  standard  PDAF  [1],  except  that  the  augmented  state  is 
conditioned  on  O^’’’  with  data  fusion  from  sensors  1  and  2.  Define  the  combined  mode-conditioned 
innovations 

mi  m2 

(28) 

0=0  6=0 
0+67^0 

Therefore,  partial  update  of  the  state  estimate 


■■=  B  {xi|9J'“,  Ml,  Jf  }  = 
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where  Kalman  gains,  axe  computed  as 

f  =  0,  for  a  =  0,  fe  =  0 

=  a  #0,6  =  0 

wt"’'  =  Pi\k-i[^  for  a  =  0,6#  0 

,Wt’'  =  Pi\k-iHt[Si]-\  for  a#0,  6#0, 
where  Hi  =  H^f'  j .  Therefore,  the  mode-conditioned  update  of  the  state  estimate 


=  < 


(30) 


mi  m2 


(31) 


0=0  6=0 
a+fe^O 


and  the  covaxiance  of 

^fclfc  ==  H\k-\  -22  -[-22 

'  '  a=0  6=0  a=0  6=0 

0-1-651^0  0+651^0 


mi  m2  .  .  ■  t  .  . 

a=0  6=0 

L  o-|-6^0 


2  2 

a=0  6=0 

J  L  a-)-65i^O 


(32) 


Step  5.  Update  of  mode  probabilities  (\/m  €  : 

4:=P[Afl|Z‘]  =  i4-Ai 

where  C  is  a  normalization  constant  such  that  EA^i  =  1- 

j 

Step  6.  Combination  of  the  mode- conditioned  estimates  (im  G  M”^  : 
The  final  augmented  state  estimate  update  at  time  k  is  given  by 

^k\k  =  Ej4|X 

and  its  covariance  is  given  by 

Pk\k  =  Xlj  l-^fclfc  [^1*  -  [^k\k  ~ 


(33) 


(34) 


(35) 


4  Simulation  Example 

The  following  example  of  tracking  a  highly  maneuvering  target  in  clutter  is  considered. 

The  True  Trajectory:  The  target  starts  at  location  [21689  10840  40]  in  Cartesian  coordinates 
in  meters.  The  initial  velocity  (in  m/s)  is  [-8.3  -399.9  0]  and  the  target  stays  at  constant  altitude 
with  a  constant  speed  of  400  m/s.  Its  trajectory  is: 
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•  a  straight  line  with  constant  velocity  between  0  and  20s, 

•  a  coordinated  tinrn  (0.15  rad/s)  with  constant  acceleration  of  60  m/s^  between  20  and  35s, 

•  a  straight  line  with  constant  velocity  between  35  and  55s, 

•  a  coordinated  turn  (0.1  rad/s)  with  constant  acceleration  of  40  m/s^  between  55  and  70s, 

•  a  straight  line  with  constant  velocity  between  70  and  90s. 

The  Target  Motion  Models:  These  are  exactly  as  in  [11].  In  each  mode  the  target  dynamics 
are  modelled  in  Cartesian  coordinates  as  r*  =  Fx/c-i  +  Gvl_^  where  the  state  of  the  target  is 
position,  velocity,  acceleration,  and  process  noise  in  each  of  the  3  Cartesian  coordinates  (x,  y,  and 
z).  Model  1  for  Nearly  constant  velocity  model  with  zero  mean  perturbation  in  acceleration;  Model 
2  for  Wiener  process  acceleration  (nearly  constant  acceleration  motion);  Model  3  for  Wiener  process 
acceleration  (model  with  large  acceleration  increments,  for  the  onset  and  termination  of  maneuvers). 
The  details  regarding  these  models  may  be  found  in  [11].  The  initial  model  probabilities  are 
g,Q  =  0.8,  /iQ  =  0.1  and  =  O-l-  The  mode  switching  probability  matrix  for  two  targets  is  also 
identical  and  is  as  in  [1]  Fig.  1  shows  the  true  trajectory  of  the  target. 

The  Sensors:  Two  sensors  (we  assume  time  synchronization  of  observations,  etc.,  but  not 
collocation)  are  used  to  obtain  the  measurements.  The  measurements  from  sensor  I  for  model  j  are 

~  +  'iiPjf  for  Z  =  1  and  2,  reflecting  range  and  azimuth  angle  for  sensor  1  (radar)  and  az¬ 

imuth  and  elevation  angles  for  sensor  2  (infrared).  The  range,  azimuth,  and  elevation  angle  transfor¬ 
mations  would  be  given  by  r/  =  {(r  -  xif  +  {y  -  gif'  -{■  {z  -  ai  =  tan~^[{y  -  yi)/{x  -  xj)], 

6/  =  tan~^[{z  —  zi)/{{x  —  x/)^  +  (y  —  yif'}^^'^],  respectively,  if  the  sensor  I  were  located  at  [x/  yi 
zi].  In  our  simulations,  we  placed  the  radar  at  [-4000  4000  0]m  and  the  infrared  sensor  at  [5000 
0  0]m.  The  measurement  noise  for  sensor  I  is  assumed  to  be  zero-mean  white  Gaussian  with 
known  covariances,  —  diag[qr,gai]  =  dmy[400m^,49mrad^]  with  qr  and  qa\  denoting  the  vari¬ 
ances  for  the  radar  range  and  azimuth  measurement  noises,  respectively,  and  =  diag[qa2,qe]  = 
diag[4mra(f,4mrad^]  with  qa2  and  qg  denoting  the  variances  for  the  infrared  sensor  azimuth  and 
elevation  measurement  noises,  respectively.  The  sampling  interval  was  T=ls  and  it  was  assumed 
that  the  probability  of  detection  Pd=0.997  for  both  sensors. 

The  Clutter:  For  generating  false  measurements  in  simulations,  the  clutter  was  assumed 
to  be  Poisson  distributed  with  expected  number  of  Ai  =  50  x  10“®/m  mrad  for  sensor  1  and 
A2  =  3.5  X  10“‘*/m^  mrad  for  sensor  2.  These  statistics  were  used  for  generating  the  clutter  in  all 
simulations.  However,  a  nonparametric  clutter  model  was  used  for  implementing  all  the  algorithms 
for  target  tracking. 

Other  Parameters:  The  gates  for  setting  up  the  validation  regions  for  both  the  sensors  were 
based  on  the  threshold  7=16  corresponding  to  a  gate  probability  Pg=0.9997  . 
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Figure  1:  True  trajectory  of  maneuvering  target  (read  left  to  right,  top  to  bottom):  (a)  Position  in 
xy  plane;  (b)  x  and  y  velocities;  (c)  x  and  y  accelerations;  (d)  magnitude  of  acceleration. 


Simulation  Results:  The  results  were  obtained  from  100  Monte  Carlo  runs.  Fig.  2  shows 
the  results  of  proposed  IMM/MSPDAF  based  on  100  runs  comparing  with  the  standard  sequential 
IMM/MSPDAF  in  terms  of  the  RMSE  (root  mean-square  error)  in  position.  It  is  seen  that  parallel 
updating  can  significantly  improve  the  accuracy  of  track  estimation.  In  addition,  it  is  also  seen 
that  tracking  with  separated  sensors  can  improve  the  accuracy  of  track  estimation  compared  to 
using  collocated  sensors. 

5  Conclusions 

We  investigated  an  IMM/MSPDAF  algorithm  with  parallel  detection  fusion  for  tracking  a  highly 
maneuvering  target  in  clutter.  The  proposed  algorithm  was  illustrated  via  a  simulation  example 
where  it  outperformed  the  IMM/MSPDAF  algorithm  with  sequential  updating  [1]. 
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Sequential  Updating  (w/  sensor  collocation)  Sequential  Updating  (w/o  sensor  collocation) 


Figure  2:  Performance  comparison  between  the  proposed  parallel  detection  fusion  IMM /MSPDAF 
and  sequential  updating  IMM/MSPDAF  of  [1]  (read  left  to  right,  top  to  bottom);  (a)  sequential 
updating  [1]  with  collocated  sensors,  (b)  sequential  updating  [1]  with  separated  sensors,  (c)  pro¬ 
posed  parallel  detection  fusion  with  collocated  sensors,  (d)  proposed  parallel  detection  fusion  with 
separated  sensors. 
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